
rm(list=ls())

setwd("workingdirectory")

library(foreign)
knutsen <- read.dta("Knutsen_Wig_CPS_final_forR.dta")

####THIS REPLICATES FIGURE 3 (to Replicate figure A2 (with expansive real-time democracy coding, just swap "dem_realtimeDD2" with "dem_realtimeDD"))


########################################################################################################################
#######                 SIMULATING THE PREDICTED ONSET OF CIVIL WAR GIVEN DIFFERENT INFORMATION SETS ###################
########################################################################################################################


###FIGURE 3 in PAPER

#######################REAL-TIME DD 

##########
#ORIGINAL#
##########
library(MASS)


civpred <- glm(acdonset2~ dem_realtimeDD2 + lgdpcapl + lpopl + gwthover5yrs + ethnfrac +  acdpeaceyears + acdspline1 + acdspline2 + acdspline3 + year , 
                    data=knutsen[knutsen$year<1970,], family=binomial(link="logit"))


beta <- coef(civpred)
covvar <- vcov(civpred)

m <- with(knutsen, cbind(1,seq(from=0, to=1, by=1),median(
  lgdpcapl,na.rm=T), median(lpopl, na.rm=T) ,median(gwthover5yrs, na.rm=T), median(ethnfrac,na.rm=T), median(acdpeaceyears,na.rm=T), median(acdspline1,na.rm=T),
  median(acdspline2,na.rm=T), median(acdspline3,na.rm=T) ,
  median(year, na.rm=T)                   
))

beta.sim <- mvrnorm(10000, beta,covvar)
xb <- m %*% t(beta.sim)
prob <- 1/(1+exp(-xb))
predprobs <- t(apply(prob,1,quantile, probs = c(.10,.5,.90)))
meanprob <- (predprobs[2,])

################
#Iteration #####
################

years <- as.vector(seq(from=1971, to=2003,by=1))

for(i in years){
civpred <- glm(acdonset2~ dem_realtimeDD2 + lgdpcapl + lpopl + gwthover5yrs + ethnfrac +  acdpeaceyears + acdspline1 + acdspline2 + acdspline3 + year , 
               data=knutsen[knutsen$year<i,], family=binomial(link="logit"))
beta <- coef(civpred)
covvar <- vcov(civpred)

m <- with(knutsen, cbind(1,seq(from=0, to=1, by=1),median(
  lgdpcapl,na.rm=T), median(lpopl, na.rm=T) ,median(gwthover5yrs, na.rm=T), median(ethnfrac,na.rm=T), median(acdpeaceyears,na.rm=T), median(acdspline1,na.rm=T),
  median(acdspline2,na.rm=T), median(acdspline3,na.rm=T) ,
  i                   
))

beta.sim <- mvrnorm(10000, beta,covvar)
xb <- m %*% t(beta.sim)
prob <- 1/(1+exp(-xb))
predprobs <- t(apply(prob,1,quantile, probs = c(.10,.5,.90)))

meanprob2 <- (predprobs[2,])
#colnames(meanprob2) <- c("predicted prob")
meanprob <- rbind(meanprob, meanprob2)
}

#######################BASELINE DD

##########
#ORIGINAL#
##########


civpred <- glm(acdonset2~ democracy + lgdpcapl + lpopl + gwthover5yrs + ethnfrac +  acdpeaceyears + acdspline1 + acdspline2 + acdspline3 + year , 
               data=knutsen[knutsen$year<1970,], family=binomial(link="logit"))

beta <- coef(civpred)
covvar <- vcov(civpred)

m <- with(knutsen, cbind(1,seq(from=0, to=1, by=1),median(
  lgdpcapl,na.rm=T), median(lpopl, na.rm=T) ,median(gwthover5yrs, na.rm=T), median(ethnfrac,na.rm=T), median(acdpeaceyears,na.rm=T), median(acdspline1,na.rm=T),
  median(acdspline2,na.rm=T), median(acdspline3,na.rm=T) ,
  median(year, na.rm=T)                   
))

beta.sim <- mvrnorm(10000, beta,covvar)
xb <- m %*% t(beta.sim)
prob <- 1/(1+exp(-xb))
predprobs <- t(apply(prob,1,quantile, probs = c(.10,.5,.90)))
meanprob3 <- (predprobs[2,])


################
#Iteration #####
################

years <- as.vector(seq(from=1971, to=2003,by=1))

for(i in years){
  civpred <- glm(acdonset2~ democracy + lgdpcapl + lpopl + gwthover5yrs + ethnfrac +  acdpeaceyears + acdspline1 + acdspline2 + acdspline3 + year , 
                 data=knutsen[knutsen$year<i,], family=binomial(link="logit"))
  beta <- coef(civpred)
  covvar <- vcov(civpred)
  
  m <- with(knutsen, cbind(1,seq(from=0, to=1, by=1),median(
    lgdpcapl,na.rm=T), median(lpopl, na.rm=T) ,median(gwthover5yrs, na.rm=T), median(ethnfrac,na.rm=T), median(acdpeaceyears,na.rm=T), median(acdspline1,na.rm=T),
    median(acdspline2,na.rm=T), median(acdspline3,na.rm=T) ,
    i                   
  ))
  
  beta.sim <- mvrnorm(10000, beta,covvar)
  xb <- m %*% t(beta.sim)
  prob <- 1/(1+exp(-xb))
  predprobs <- t(apply(prob,1,quantile, probs = c(.10,.5,.90)))
  
  meanprob4 <- (predprobs[2,])
  meanprob3 <- rbind(meanprob3, meanprob4)
}

years <- as.data.frame(years)
years <- rbind(1970, years)
preds <- as.data.frame(cbind(years, meanprob, meanprob3))

colnames(preds) <- c("year", "c1r", "realtimeDD","c2r", "c1b", "baselineDD", "c2b")


pdf("civwarpreds2.pdf")
plot(1,1,ylim=c(0,0.10),xlim=c(1970,2003), type="n",
     xlab="year", ylab="Predicted probability of civil war onset in a democracy")
library(denstrip)
polygon(c(preds$year
          , rev(preds$year))
        , c(preds$c1r
            , rev(preds$c2r))
        , col = 'grey'
        , border = NA)
lines(1970:2003, preds$realtimeDD, col=2)
lines(1970:2003, preds$baselineDD, col=1)
legend("bottomright", legend=c("DD (real time)", "DD"),col=c(2,1), lty=1)
dev.off()


########################################################################################################################
#######                 SIMULATING THE PREDICTED ONSET OF COUP ATTEMPTS GIVEN DIFFERENT INFORMATION SETS ###################
########################################################################################################################

#######################REAL-TIME DD 

##########
#ORIGINAL#
##########
library(MASS)

coupred <- glm(attempt ~ dem_realtimeDD2 +  milreg + chmilex + lmilper +chgdp + lgdppcl + instab 
                    + powthy_peace + spline1 + spline2 + spline3
                    , data=knutsen[knutsen$year<1980,], family=binomial(link="logit"))


beta <- coef(coupred)
covvar <- vcov(coupred)

m <- with(knutsen, cbind(1,seq(from=0, to=1, by=1),median(milreg, na.rm=T),
  median(chmilex, na.rm=T), median(lmilper, na.rm=T), median(chgdp, na.rm=T), median(lgdppcl, na.rm=T), median(instab, na.rm=T), 
         median(powthy_peace, na.rm=T), median(spline1, na.rm=T), median(spline2, na.rm=T) , median(spline3, na.rm=T)                 
))

beta.sim <- mvrnorm(10000, beta,covvar)
xb <- m %*% t(beta.sim)
prob <- 1/(1+exp(-xb))
predprobs <- t(apply(prob,1,quantile, probs = c(.10,.5,.90)))
meanprob <- (predprobs[2,])

################
#Iteration #####
################

years <- as.vector(seq(from=1981, to=2003,by=1))

for(i in years){
  coupred <- glm(attempt ~ dem_realtimeDD2 +  milreg + chmilex + lmilper +chgdp + lgdppcl + instab 
                 + powthy_peace + spline1 + spline2 + spline3
                 , data=knutsen[knutsen$year<i,], family=binomial(link="logit"))
  beta <- coef(coupred)
  covvar <- vcov(coupred)
  
  m <- with(knutsen, cbind(1,seq(from=0, to=1, by=1),median(milreg, na.rm=T),
                           median(chmilex, na.rm=T), median(lmilper, na.rm=T), median(chgdp, na.rm=T), median(lgdppcl, na.rm=T), median(instab, na.rm=T), 
                           median(powthy_peace, na.rm=T), median(spline1, na.rm=T), median(spline2, na.rm=T) , median(spline3, na.rm=T)                 
  ))
  
  
  beta.sim <- mvrnorm(10000, beta,covvar)
  xb <- m %*% t(beta.sim)
  prob <- 1/(1+exp(-xb))
  predprobs <- t(apply(prob,1,quantile, probs = c(.10,.5,.90)))
  
  meanprob2 <- (predprobs[2,])
  #colnames(meanprob2) <- c("predicted prob")
  meanprob <- rbind(meanprob, meanprob2)
}

#######################BASELINE DD

##########
#ORIGINAL#
##########


coupred <- glm(attempt ~ democracy +  milreg + chmilex + lmilper +chgdp + lgdppcl + instab 
               + powthy_peace + spline1 + spline2 + spline3
               , data=knutsen[knutsen$year<1980,], family=binomial(link="logit"))


beta <- coef(coupred)
covvar <- vcov(coupred)

m <- with(knutsen, cbind(1,seq(from=0, to=1, by=1),median(milreg, na.rm=T),
                         median(chmilex, na.rm=T), median(lmilper, na.rm=T), median(chgdp, na.rm=T), median(lgdppcl, na.rm=T), median(instab, na.rm=T), 
                         median(powthy_peace, na.rm=T), median(spline1, na.rm=T), median(spline2, na.rm=T) , median(spline3, na.rm=T)                 
))
beta.sim <- mvrnorm(10000, beta,covvar)
xb <- m %*% t(beta.sim)
prob <- 1/(1+exp(-xb))
predprobs <- t(apply(prob,1,quantile, probs = c(.10,.5,.90)))
meanprob3 <- (predprobs[2,])


################
#Iteration #####
################

for(i in years){
  coupred <- glm(attempt ~ democracy +  milreg + chmilex + lmilper +chgdp + lgdppcl + instab 
                 + powthy_peace + spline1 + spline2 + spline3
                 , data=knutsen[knutsen$year<i,], family=binomial(link="logit"))
  beta <- coef(coupred)
  covvar <- vcov(coupred)
  
  m <- with(knutsen, cbind(1,seq(from=0, to=1, by=1),median(milreg, na.rm=T),
                           median(chmilex, na.rm=T), median(lmilper, na.rm=T), median(chgdp, na.rm=T), median(lgdppcl, na.rm=T), median(instab, na.rm=T), 
                           median(powthy_peace, na.rm=T), median(spline1, na.rm=T), median(spline2, na.rm=T) , median(spline3, na.rm=T)                 
  ))
  
  
  beta.sim <- mvrnorm(10000, beta,covvar)
  xb <- m %*% t(beta.sim)
  prob <- 1/(1+exp(-xb))
  predprobs <- t(apply(prob,1,quantile, probs = c(.10,.5,.90)))
  
  meanprob4 <- (predprobs[2,])
  meanprob3 <- rbind(meanprob3, meanprob4)
}

years <- as.data.frame(years)
years <- rbind(1980, years)
preds <- as.data.frame(cbind(years, meanprob, meanprob3))

colnames(preds) <- c("year", "c1r", "realtimeDD","c2r", "c1b", "baselineDD", "c2b")


pdf("couppreds2.pdf")
plot(1,1,ylim=c(0,0.03),xlim=c(1980,2003), type="n",
     xlab="year", ylab="Predicted probability of a coup attempt in a democracy")
library(denstrip)
polygon(c(preds$year
          , rev(preds$year))
        , c(preds$c1r
            , rev(preds$c2r))
        , col = 'grey'
        , border = NA)
lines(1980:2003, preds$realtimeDD, col=2)
lines(1980:2003, preds$baselineDD, col=1)
legend("bottomright", legend=c("DD (real time)", "DD"),col=c(2,1), lty=1)
dev.off()
